function X = myim2col(A, varargin)
% myim2col        Rearrange image blocks into columns.
%                 This is more flexible variant of Matlab im2col 
% The most important feature is that the columns can be made from a
% coefficient domain of the image, i.e. a transform is done first then the
% blocks are extracted.
% The 'inverse' is mycol2im.m
% The implemetation use some few functions besides standard matlab
% functions, i.e: myimadjust, mylwt2, getLOT and getELT
%
% X = myim2col(A, argName, argVal, ...);
%----------------------------------------------------------------------
%  A     an image, ex: A = double(imread('lena.bmp')) - 128;
%  There may be an additional number of input arguments, a struct, a cell 
%        or as pairs: argName, argVal, ...
%  'o' or 'offset' which may be [0,0] (default) or other values [Mo,No]
%        This is done first, i.e. the first operation on the image is to
%        remove the Mo first rows and the No first columns, thus making
%        A(Mo+1,No+1) as the new (upper left) corner of the image
%        When used, offset reduce the size of the image
%  't' or 'transform'. 
%        The option value may be the name of a transform: 'dct', 'lot' 
%        or 'elt', or the transform given as a matrix. Note that in
%        analysis part (here) the transposed of this matrix is used, while
%        for synthesis (mycol2im) the matrix is used as given, this means
%        that if the transform is orthogonal the same matrix can be given
%        both in myim2col and mycol2im. Alternatively, the option value
%        can be the name of a wavelet or the wavelet lifting scheme given 
%        directly as a cell array. Wavelt name as in mylwt2, that is as
%        liftwave (in Matlab  Wavelet Toolbox) or 'j97' or 'm97'.
%        The transform is always a separable transform, and except for
%        'dct', 'lot' and 'elt' with same size in both directions.
%        Default option value is 'none' for reshape in pixel domain.
%  's' or 'size' is for the size of the transform. Only for 'dct', 'lot'
%        and 'elt' it may be different in the two dimension, [Ms,Ns].
%        For wavelet Ns = Ms, and it should be 2,4,8,16 or 32.
%  'a' or 'adjust'. Before the transform is done the image size must be
%        adjusted to match the size of the transform. This is done by
%        myimadjust.m and the value of the adjust parameter is used as
%        method in myimadjust: 'none' or 'extend' is preferred.
%  'n' or 'neighborhood' can be used to extract special neighborhoods
%        around a pixel. The option value may be 8, [8,8], ones(8) or a 
%        (logical)  matrix indicating which elements to include, 
%        for example: kron( ones(4), [1,0]'*[1,0] );
%        default is size of transform or [4,4] if transform is 'none'.
%  'i' or 'increment' which may be  'distinct' for (size of) neighborhood
%        blocks (default), or 'sliding' for [1,1] or numbers [Mi,Ni]
%  'v' or 'verbose' to indicate verboseness
%----------------------------------------------------------------------
% examples:  (see test01 for mere examples)
%  A = double(imread('lena.bmp')) - 128;
%  param = struct('t','m97', 's',[8,8], 'a','extend', 'i','distinct', 'v',1);
%  X = myim2col(A, param);
% to do just the transform (in-place) use for example:
%  X = myim2col(A, 't','m97', 's',8, 'n',[size(A,1),1], 'v',1);

%----------------------------------------------------------------------
% Copyright (c) 2009.  Karl Skretting.  All rights reserved.
% University of Stavanger (Stavanger University), Signal Processing Group
% Mail:  karl.skretting@uis.no   Homepage:  http://www.ux.his.no/~karlsk/
% 
% HISTORY:  dd.mm.yyyy
% Ver. 1.0  26.11.2009  Made function
% Ver. 1.1  08.03.2010  KS: made the function 'cleaner'
% Ver. 1.2  11.05.2011  KS: corrected an error for I and J when blocks are found  
%----------------------------------------------------------------------

%% check if imwrite should be called 
if (nargin < 1)
    error('myim2col: too few input arguments.');
end

%% default options
Mo=0; No=0;   % offset
tr = 'none';  % transform
Ms=0; Ns=0;   % size (of transform)
amet = '';    % method in myimadjust
nei = [];     % neighborhood, input or given by tr
Mn=0; Nn=0;   % size of nei (neighborhood)
Mi=0; Ni=0;   % increment
verbose = 0;
lotrho = 0.95;  % these three must be the same in mycol2im
eltrho = 0.95;
eltarg2 = 0.7;

%% get the options
nofOptions = nargin-1;
optionNumber = 1;
fieldNumber = 1;
while (optionNumber <= nofOptions)
    if isstruct(varargin{optionNumber})
        sOptions = varargin{optionNumber}; 
        sNames = fieldnames(sOptions);
        opName = sNames{fieldNumber};
        opVal = sOptions.(opName);
        % next option is next field or next (pair of) arguments
        fieldNumber = fieldNumber + 1;  % next field
        if (fieldNumber > numel(sNames)) 
            fieldNumber = 1;
            optionNumber = optionNumber + 1;  % next pair of options
        end
    elseif iscell(varargin{optionNumber})
        sOptions = varargin{optionNumber}; 
        opName = sOptions{fieldNumber};
        opVal = sOptions{fieldNumber+1};
        % next option is next pair in cell or next (pair of) arguments
        fieldNumber = fieldNumber + 2;  % next pair in cell
        if (fieldNumber > numel(sOptions)) 
            fieldNumber = 1;
            optionNumber = optionNumber + 1;  % next pair of options
        end
    else
        opName = varargin{optionNumber};
        opVal = varargin{optionNumber+1};
        optionNumber = optionNumber + 2;  % next pair of options
    end
    % interpret opName and opVal
    if (strcmpi(opName,'offset') || strcmpi(opName,'o'))
        if ((numel(opVal) == 1) && isnumeric(opVal))
            Mo = max(0, floor(opVal));
            No = Mo;
        elseif ((numel(opVal) == 2) && isnumeric(opVal))
            Mo = max(0, floor(opVal(1)));
            No = max(0, floor(opVal(2)));
        else
            disp('myim2col: illegal option offset. We ignore it.');
        end
    end
    if (strcmpi(opName,'transform') || strcmpi(opName,'t'))
        tr = opVal;
    end
    if (strcmpi(opName,'size') || strcmpi(opName,'s'))
        if ((numel(opVal) == 1) && isnumeric(opVal))
            Ms = max(0, floor(opVal));
            Ns = Ms;
        elseif ((numel(opVal) == 2) && isnumeric(opVal))
            Ms = max(0, floor(opVal(1)));
            Ns = max(0, floor(opVal(2)));
        else
            error('myim2col: illegal option size. We ignore it.');
        end
    end
    if (strcmpi(opName,'adjust') || strcmpi(opName,'a'))
        if (strcmpi(opVal,'none') || strcmpi(opVal,'extend') || ...
            strcmpi(opVal,'zeros') || strcmpi(opVal,'periodic') || ...    
            strcmpi(opVal,'mirror') )
            amet = opVal;
        else
            disp('myim2col: illegal option for adjust method. We ignore it.');
        end
    end
    if (strcmpi(opName,'neighborhood') || strcmpi(opName,'n'))
        if ((numel(opVal) <= 2) && isnumeric(opVal))
            nei = ones(opVal);
            [Mn, Nn] = size(nei);
        elseif ((numel(opVal) > 2) && isnumeric(opVal))
            nei = opVal;
            [Mn, Nn] = size(nei);
        else
            error('myim2col: illegal option neighborhood. We ignore it.');
        end
    end
    if (strcmpi(opName,'increment') || strcmpi(opName,'i'))
        if ((numel(opVal) == 1) && isnumeric(opVal))
            Mi = max(0, floor(opVal));
            Ni = Mi;
        elseif ((numel(opVal) == 2) && isnumeric(opVal))
            Mi = max(0, floor(opVal(1)));
            Ni = max(0, floor(opVal(2)));
        elseif strcmpi(opVal, 'distinct');
            Mi = 0;
            Ni = 0;
        elseif strcmpi(opVal, 'sliding');
            Mi = 1;
            Ni = 1;
        else
            error('myim2col: illegal option increment. We ignore it.');
        end
    end
    if strcmpi(opName,'verbose') || strcmpi(opName,'v')
        if (islogical(opVal) && opVal); verbose = 1; end;
        if isnumeric(opVal); verbose = opVal(1); end;
    end
end

%% arguments are now read and interpreted. Start processing
% if a special data set is wanted, the rest of the arguments are ignored


%% check and display options
if ((Ms == 0) || (Ns == 0))
    if strcmpi(tr,'none')
        Ms = 4; Ns = 4;
    else
        Ms = 8; Ns = 8;
    end
end
if (~iscell(tr) && ((strcmpi(tr,'lot') || strcmpi(tr,'elt'))))
    Ms = Ms - mod(Ms,2);   % Ms and Ns is made even
    Ns = Ns - mod(Ns,2);   
end
if ( iscell(tr) || ... 
        ( ischar(tr) && ... 
          ~strcmpi(tr,'none') && ...
          ~strcmpi(tr,'lot') && ...
          ~strcmpi(tr,'elt') && ...
          ~strcmpi(tr,'dct') )  )
    % a wavelet as in mylwt2
    Ms = max(2,2^floor(log2(Ms)));
    Ns = Ms;
end
if (isnumeric(tr)) % transform is given as a matrix
    Ms = size(tr,2);
    Ns = Ms;
    P = floor(size(tr,1)/Ms);  % overlap factor, should be integer
    if (Ms*P ~= size(tr,1))
        tr = tr(1:(Ms*P));
    end
end

if (numel(amet) == 0)
    if strcmpi(tr,'none')
        amet = 'none';
    else
        amet = 'extend';
    end
end
if ((Mn == 0) || (Nn == 0))
    nei = ones(Ms,Ns);
    [Mn, Nn] = size(nei);
end
if ((Mi == 0) || (Ni == 0))
    Mi = Mn; 
    Ni = Nn;
end
if verbose
    disp(['myim2col: Input image is ',int2str(size(A,1)),...
        'x',int2str(size(A,2))]);
end
        

%% remove what is left and above offset, and perhaps adjust the size
if ((Mo > 0) || (No > 0))
    if verbose
        disp(['Remove ',int2str(Mo),' left columns and ',...
            int2str(No),' top rows from image.']);
    end
    A = A((Mo+1):end, (No+1):end);
end
if ~strcmpi(tr,'none')  
    if verbose
        disp(['Adjust image to fit ',int2str(Ms),'x',...
            int2str(Ns),' blocks by method ',amet,'.']);
    end
    A = myimadjust(A, amet, [Ms, Ns]);
end
[M,N] = size(A);  % the new adjusted image size
if verbose
    disp(['Adjusted image is ',int2str(M),...
        'x',int2str(N)]);
end

%% do the transform
if strcmpi(tr,'none')
    % do nothing
elseif strcmpi(tr,'dct')
    if verbose
        disp(['Do DCT with size ',int2str(Ms),'x',int2str(Ns),'.']);
    end
    if (Ms > 1)
        A = dct( reshape(A,Ms,M*N/Ms) );    % the columns (of each block)
        A = reshape(A,M,N);  % back to size MxN
    end
    if (Ns > 1)
        A = dct( reshape(A',Ns,M*N/Ns) ); % the rows
        A = reshape(A,N,M)';  % back to size MxN
    end
elseif strcmpi(tr,'lot')
    if verbose
        disp(['Do LOT with size ',int2str(Ms),'x',int2str(Ns),'.']);
    end
    if (Ms > 1)
        F = getLOT(Ms, lotrho);  % get the synthesis vectors of LOT
        F1 = F(1:Ms,:);
        F2 = F((Ms+1):(2*Ms),:);
        A = F1' * reshape(A,Ms,M*N/Ms) + ...
            F2' * reshape([A((Ms+1):M,:);A(1:Ms,:)],Ms,M*N/Ms);
        A = reshape(A,M,N);
    end
    if (Ns > 1)
        F = getLOT(Ns, lotrho);
        F1 = F(1:Ns,:);
        F2 = F((Ns+1):(2*Ns),:);
        A = F1' * reshape(A',Ns,M*N/Ns) + ...
            F2' * reshape([A(:,(Ns+1):N),A(:,1:Ns)]',Ns,M*N/Ns);
        A = reshape(A,N,M)';
    end
elseif strcmpi(tr,'elt')
    if verbose
        disp(['Do ELT with size ',int2str(Ms),'x',int2str(Ns),'.']);
    end
    if (Ms > 1)
        F = reshape(getELT(Ms,eltarg2,eltrho)',Ms,Ms,4); 
        A = F(:,:,1) * reshape(A,Ms,M*N/Ms) + ...
            F(:,:,2) * reshape([A((Ms+1):M,:);A(1:Ms,:)],Ms,M*N/Ms) + ...
            F(:,:,3) * reshape([A((2*Ms+1):M,:);A(1:(2*Ms),:)],Ms,M*N/Ms) + ...
            F(:,:,4) * reshape([A((3*Ms+1):M,:);A(1:(3*Ms),:)],Ms,M*N/Ms);
        A = reshape(A,M,N);
    end
    if (Ns > 1)
        F = reshape(getELT(Ns,eltarg2,eltrho)',Ns,Ns,4); 
        A = F(:,:,1) * reshape(A',Ns,M*N/Ns) + ...
            F(:,:,2) * reshape([A(:,(Ns+1):N),A(:,1:Ns)]',Ns,M*N/Ns) + ...
            F(:,:,3) * reshape([A(:,(2*Ns+1):N),A(:,1:(2*Ns))]',Ns,M*N/Ns) + ...
            F(:,:,4) * reshape([A(:,(3*Ns+1):N),A(:,1:(3*Ns))]',Ns,M*N/Ns);
        A = reshape(A,N,M)';
    end
elseif isnumeric(tr)   % tr is (Ms*P x Ns)
    if P==1  % no overlap
        F = tr';   % Ms x Ms  (= Ns x Ns), the transposed
        if verbose
            disp(['Do supplied transform with size ',...
                int2str(Ms),'x',int2str(Ms),'.']);
        end
    else   
        F = reshape(tr',Ms,Ms,P);  % Ms x Ms x P  % the transposed
        if verbose
            disp(['Do supplied overlapping transform with size ',...
                int2str(Ms),'x',int2str(Ms),' (P=',int2str(P),').']);
        end
    end
    % transform columns
    % disp([size(F),size(A),M,N,Ms,Ns,P]);
    Y = F(:,:,1) * reshape(A,Ms,M*N/Ms);
    for p=2:P
        Y = Y + F(:,:,p) * reshape([A(((p-1)*Ms+1):M,:);A(1:((p-1)*Ms),:)],Ms,M*N/Ms);
    end
    A = reshape(Y,M,N);
    % transform rows
    Y = F(:,:,1) * reshape(A',Ns,M*N/Ns);
    for p=2:P
        Y = Y + F(:,:,p) * reshape([A(:,((p-1)*Ns+1):N),A(:,1:((p-1)*Ns))]',Ns,M*N/Ns);
    end
    A = reshape(Y,N,M)';
    %
elseif (iscell(tr) || ischar(tr)) % a wavelet as in mylwt2
    if verbose
        if (iscell(tr)); temp = 'by given liftwave'; 
        else temp = tr;
        end
        disp(['Do wavelet ',temp,', size ',int2str(Ms),'x',int2str(Ns),'.']);
    end
    A = mylwt2(A, tr, log2(Ms));
end

%% find the blocks
index = find(nei);
I = 1:Mi:(size(A,1)-Mn+1);   % was I = (Mo+1):Mi:... before May 10 2011, but that must be error 
J = 1:Ni:(size(A,2)-Nn+1);   % and here 
L = numel(I)*numel(J);
if verbose
    disp(['neighborhood is ',int2str(numel(index)),' pixels from ',...
        int2str(Mn),'x',int2str(Nn),' block.']);
    disp(['Get ',int2str(L),' columns using ',...
        'increment step ',int2str(Mi),' and ',int2str(Ni),'.']);
end
X = zeros(numel(index), L);
m = Mn-1; 
n = Nn-1;
k = 1;
for j=J
    for i=I
        block = A(i:(i+m),j:(j+n));
        X(:,k) = block(index);
        k = k+1;
    end
end

return


